Human Colon Xenium in situ ST Dataset, Nuclei Segmentation

First we need to set up the environment and load the packages we will use for this workshop.

library(Seurat): Loads the Seurat package, which is a comprehensive toolkit for single-cell RNA sequencing and spatial transcriptomics data analysis. It provides a wide range of functions for data preprocessing, normalization, clustering, dimensionality reduction, and visualization. Explore documentation here: https://satijalab.org/seurat/

library(ggplot2): Loads the ggplot2 package, a powerful and flexible system for creating static visualizations in R. Explore documentation here: https://ggplot2.tidyverse.org/

library(scCustomize): Loads the scCustomize package, which provides custom functions and themes to enhance the visualization and analysis capabilities of single-cell and spatial transcriptomics data, often in conjunction with Seurat. Explore documentation here: https://samuel-marsh.github.io/scCustomize/

library(readr): Loads readr package for fast and friendly reading of rectangular data, such as CSV files, into R.

library(pheatmap): Loads pheatmap package, which is for creating pretty heatmaps, offering better control over heatmap customization compared to base R.

library(matrixStats): matrixStats provides highly optimized functions for matrix operations, particularly useful for computing row and column summaries.

library(spdep): spdep stands for Spatial Dependence and Spatial Autocorrelation, and it provides functions for spatial data analysis, including spatial weights generation, spatial autocorrelation statistics, and spatial regression.

library(geojsonR) The geojsonR library is used for handling GeoJSON data in R. GeoJSON is a format for encoding a variety of geographic data structures using JavaScript Object Notation (JSON). It is sometimes used as a format for storing cell segmentation boundaries.

library(Seurat)

Loading required package: SeuratObject
Loading required package: sp
The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
which was just loaded, will retire in October 2023.
Please refer to R-spatial evolution reports for details, especially
https://r-spatial.org/r/2023/05/15/evolution4.html.
It may be desirable to make the sf package available;
package maintainers should consider adding sf to Suggests:.
The sp package is now running under evolution status 2
     (status 2 uses the sf package in place of rgdal)

Attaching package: ‘SeuratObject’

The following objects are masked from ‘package:base’:

    intersect, t

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

library(ggplot2)
library(scCustomize)

scCustomize v2.1.2
If you find the scCustomize useful please cite.
See 'samuel-marsh.github.io/scCustomize/articles/FAQ.html' for citation info.

library(readr)
library(pheatmap)
library(matrixStats)
library(spdep)

Loading required package: spData
To access larger datasets in this package, install the spDataLarge
package with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`
Loading required package: sf
Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE

library(geojsonR)

Sets the path to the directory containing the Xenium output data - this is the directory where all of the outputs are stored.

data_dir <- "/project/shared/spatial_data_camp/datasets/DATASET2/XENIUM_COLORECTAL_CANCER"

ReadXenium reads Xenium spatial transcriptomics data from a specified directory using a Seurat wrapper function that supports this data format. Xenium data typically includes expression matrices and spatial coordinates, along with other information about cell centroids and segmentations and coordinates of individual transcripts.

data_dir: The path to the directory containing the Xenium data, set in the previous step. outs = c(“matrix”, “microns”): Specifies the outputs to read from the data directory. matrix refers to summarised cell by gene matrix and microns refers to individual transcript coordinates.

type = c(“centroids”, “segmentations”): Indicates the types of spatial information to include - here, we are reading ib both cell centroid coordinates and cell boundary segmentations.

data <- ReadXenium(data_dir, outs = c("matrix", "microns"), type=c("centroids", "segmentations"))

10X data contains more than one type and is being returned as a list containing matrices of each type.
|--------------------------------------------------|
|==================================================|

This provides us a list of data:

names(data)

[1] "matrix"        "microns"       "centroids"     "segmentations"

Matrix is further split into gene expression matrix and various control probes and codewords. Different platforms and platform versions include different control probes. As this will vary, it’s important to check and understand what the specific controls in your own data are.

Here, negative control probes are probes that are added to the reaction but target non-biological sequences and should not bind any tissue RNA. Negative control codewords are valid codewords, but no probes with that codeword added to the reaction. This effectively tells us how good the transcript calling algorithm is.

names(data$matrix)

[1] "Gene Expression"           "Negative Control Codeword"
[3] "Negative Control Probe"    "Unassigned Codeword"      

Read in additional information about the cells - this gives us pre-calculated information, for example segmented cell or nucleus size for each cell.

cell_meta_data <- read.csv(file.path(data_dir, "cells.csv.gz"))
rownames(cell_meta_data) <- cell_meta_data$cell_id
head(cell_meta_data)

We will start by creating a basic seurat object from the data.

CreateSeuratObject function initializes a Seurat object using the provided gene expression matrix and optional metadata.

counts: The gene expression matrix, which contains the raw count data for each gene in each cell. data$matrix[[“Gene Expression”]]: Specifies the gene expression matrix extracted from the loaded Xenium data. Here, we leave out the control probes for now.

assay: The name of the assay - you can call it anything you like. Here, we go with “XENIUM”.

meta.data: Metadata associated with the cells or spots. Here, we add the cell statistics we read in earlier as cell_meta_data.

By printing the seurat object, we can see that we read in ~ 30,000 cells with measures for 325 genes

seurat <- CreateSeuratObject(counts = data$matrix[["Gene Expression"]],

Warning message:
In png(..., res = dpi, units = "in") :
  unable to open connection to X11 display ''

                                 assay = "XENIUM",
                                 meta.data = cell_meta_data)
seurat

An object of class Seurat 
325 features across 647524 samples within 1 assay 
Active assay: XENIUM (325 features, 0 variable features)
 1 layer present: counts

Adding spatial coordinates to a Seurat object allows for spatially resolved analysis and visualization. This requires creating objects for centroids and segmentations we read in earlier, and then integrating these with the main Seurat object.

CreateFOV: This function creates a field of view (FOV) object that includes spatial information about the centroids, segmentations, and molecule coordinates. An FOV can be the entire slide, or a selected region within a slide - i.e. it does not need to have entries for all the cells in the seurat object.

coords: A list containing the centroids and/or segmentation data. For larger datasets, it can be quicker to only load centroids, as this minimises the amount of data points.

centroids = CreateCentroids(data\(centroids)*: Creates a centroids object from the centroid data in the Xenium dataset. *segmentation = CreateSegmentation(data\)segmentations): Creates a segmentation object from the segmentation data in the Xenium dataset.

type = c(“segmentation”, “centroids”): Specifies the types of spatial data being included, which are segmentation and centroid data.

molecules = data$microns: The spatial coordinates of individual transcripts/molecules in the data. This is optional - for larger datasets, skipping transcript coordinates can be a good idea.

seurat[[“COLON”]] <- coords: Adds the created FOV object to the Seurat object under the new FOV name “COLON”. This can be named (almost) anything - but, avoid using underscores as this can create some unexpected behaviours later.

TIP: LoadXenium() is a wrapper that would load in both cell counts matrix and spatial coordinates in one function, simplifying these steps. However, in situ platforms are evolving at a very fast rate and there are constant changes on how the data is stored, in particular for file formats for cell segmentation and coordinates. Here, we have broken down the steps to show how to assemble an in situ seurat object from the key components, in case the platform specific readers don’t work for your specific data.

coords <- CreateFOV(coords = list(centroids = CreateCentroids(data$centroids), 
                                  segmentation = CreateSegmentation(data$segmentations)),
                    type = c("segmentation", "centroids"),
                    assay = "XENIUM")
seurat[["COLO_CAN"]] <- coords  

Inspect the object - now, you can see we have added a spatial field of view:

seurat

An object of class Seurat 
325 features across 647524 samples within 1 assay 
Active assay: XENIUM (325 features, 0 variable features)
 1 layer present: counts
 2 spatial fields of view present: COLO_CAN CRC1

#crop the image to check which coordinates are needed
cropped <- Crop(seurat[["COLO_CAN"]], x = c(4000, 2000), y = c(1000, 4000), coords = "plot")
seurat[["CRC1"]] <- cropped

Warning: Key ‘XENIUM_’ taken, using ‘crc1_’ instead

ImageDimPlot(seurat,fov = "CRC1", axes = TRUE)

seurat_CRC1 <- seurat[,seurat$x_centroid < 4000 & seurat$x_centroid > 2000 &
                        seurat$y_centroid <4000 & seurat$y_centroid > 1000]

Warning: Not validating Centroids objectsWarning: Not validating Centroids objectsWarning: Not validating FOV objectsWarning: Not validating Centroids objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating Seurat objectsWarning: Not validating Centroids objectsWarning: Not validating Centroids objectsWarning: Not validating FOV objectsWarning: Not validating Centroids objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating Seurat objects

ImageDimPlot(seurat_CRC1, axes=TRUE)

Adding control probes and codewords as separate assays in the Seurat object allows for the tracking and analysis of technical artifacts and noise within your spatial transcriptomics data, while keeping these outputs separate from the main biological gene expression values.

Unassigned codewords are unused codewords. There is no probe in a particular gene panel that will generate the codeword.

Negative control probes are probes that exist in the panels but target non-biological sequences. They can be used to assess the specificity of the assay.

Negative control codewords are codewords in the codebook that do not have any probes matching that code. They are chosen to meet the same requirements as regular codewords and can be used to assess the specificity of the decoding algorithm.

seurat_CRC1[["Negative.Control.Codeword"]] <- CreateAssayObject(counts = data$matrix[["Negative Control Codeword"]][, WhichCells(seurat_CRC1)])

Warning: Feature names cannot have underscores ('_'), replacing with dashes ('-')Warning: Feature names cannot have underscores ('_'), replacing with dashes ('-')

seurat_CRC1[["Negative.Control.Probe"]] <- CreateAssayObject(counts = data$matrix[["Negative Control Probe"]][, WhichCells(seurat_CRC1)])

Warning: Feature names cannot have underscores ('_'), replacing with dashes ('-')Warning: Feature names cannot have underscores ('_'), replacing with dashes ('-')

seurat_CRC1[["Unassigned.Codeword"]] <- CreateAssayObject(counts = data$matrix[["Unassigned Codeword"]][, WhichCells(seurat_CRC1)])

Warning: Feature names cannot have underscores ('_'), replacing with dashes ('-')Warning: Feature names cannot have underscores ('_'), replacing with dashes ('-')

Inspect the object:

seurat_CRC1

An object of class Seurat 
541 features across 14054 samples within 4 assays 
Active assay: XENIUM (325 features, 0 variable features)
 1 layer present: counts
 3 other assays present: Negative.Control.Codeword, Negative.Control.Probe, Unassigned.Codeword
 1 spatial field of view present: COLO_CAN

Let’s start with some basic QC and visualisation of the data.

In Seurat, in situ spatial transcriptomics counterpart functions to ‘SpatialDimPlot’ and ‘SpatialFeaturePlot’ we covered yesterday are called ‘ImageFeaturePlot’ and ‘ImageDimPlot’. These have additional functionality to plot cell segmentations and individual transcript coordinates, but otherwise function exactly the same as the sequencing based ST counterparts.

First, lets visualise the total transcripts detected per cell.

As in scRNA-Seq data, this is the most basic measure of overall signal and how well the data looks.

Unlike in scRNA-Seq data or unbiased sequencing-based ST, these measures are also very heavily dependent not only on the total RNA quantity of each cell and tissue quality, but also on the target panel used for the experiment. Under-represented cell types will naturally yield fewer transcripts. Finally, the quality of cell segmentation also plays a role.

In this case, we can see that there are areas with higher and lower total transcripts detected.

Understanding your tissue and target panel here is important to delineate where these differences are biological and where they may be technical.

ImageFeaturePlot(seurat_CRC1, "nCount_XENIUM") + scale_fill_viridis_c()

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

Similarly, we can visualise the total number of gene detected per cell. You can see that this is a bit less variable across tissue.

ImageFeaturePlot(seurat_CRC1, "nFeature_XENIUM") + scale_fill_viridis_c()

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

This code examines the distribution of the number of features (genes) detected per cell in the Seurat object using a density plot and calculates specific quantiles of this distribution. This is important for understanding the variability and distribution of detected features, which can help identify potential issues such as low-quality cells and determine any filtering thresholds that may need to be applied.

If you’re coming from scRNA-Seq work, these low numbers probably look very alarming. How can you possibly work with 31 median genes per cell?

Unlike scRNA-Seq data and sequencing-based ST, both gene dropouts and noise are much, much lower in in situ ST data.

We are also working with 100-fold fewer targetted genes.

ggplot(seurat[[]], aes(nFeature_XENIUM)) + geom_density()

quantile(seurat$nFeature_XENIUM, c(0.01, 0.1, 0.5, 0.9, 0.99))

 1% 10% 50% 90% 99% 
  5  14  32  54  73 

Using ImageFeaturePlot to visualize the cell area in spatial transcriptomics data allows us to examine the spatial organization and potential heterogeneity of cell sizes within your tissue sample.

Why do we get such a difference in spatial distribution of cell sizes?

This could be due to biological differences between small and large cells - e.g. small cells like T-cells.

However, here the signal correlates with areas of low cellularisation. Therefore, it is likely this is an artefact of nuclei expansion in cell segmentation.

What is Nuclei Expansion?

Nuclei expansion in cell segmentation refers to the process of enlarging the segmented nuclei regions to approximate the boundaries of the entire cells. This technique is used to better represent the actual cell boundaries when only the nuclei have been explicitly segmented/we only have DAPI and no additional cell boundary staining. The primary goal is to provide a more accurate estimation of the cellular area, which is crucial for various downstream analyses in spatial transcriptomics and single-cell studies. In this case, nuclei expansion is constrained either by maximum distance or other nearby cells - so, where there are no other nearby cells to “bump into”, the expansion generates artificially bigger cells.

ImageFeaturePlot(seurat_CRC1, "cell_area", axes = T) + scale_fill_viridis_c()

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

We can further check that this is likely the case by plotting the ratio between nuclei and total cell area. We can see that there is a very big decrease in percentage of cell area occupied by nucleus in areas of low cell density.

The cell-to-nucleus area ratio can also potentially provide insights into cell morphology, cell type and potential changes in cellular states or conditions. For example, T-Cells can often be quite well identified by this variable alone, as they have a small cytoplasm volume. However, without a cell boundary stain, this metric mainly captures segmentation artefacts, so be careful about over-interpretation!

seurat_CRC1$cell_nucleus_ratio <- seurat_CRC1$nucleus_area / seurat_CRC1$cell_area
ImageFeaturePlot(seurat_CRC1, "cell_nucleus_ratio") + scale_fill_viridis_c()

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

If we look at the distribution, we see that we have a big tail end of overly large cells.

ggplot(seurat_CRC1[[]], aes(cell_area)) + geom_density()

In this case, we can see that as expected, there is generally a correlation between cell area and transcript detection rate.

However, we also have a group of cells where this is not the case - very large cells but relatively few transcripts. These cells are mainly submucosal stromal cells which are very poorly covered by the panel 10x have used.

ggplot(seurat_CRC1[[]], aes(nCount_XENIUM, cell_area)) + geom_point() 

We can create a filter to remove the overly large cells from the analysis.

quantile(seurat$cell_area, 0.99): Calculates the 99th percentile of the cell_area values in the Seurat object. This value serves as a threshold to identify the largest 1% of cells - but what is a sensible threshold, if any, depends on your tissue.

seurat\(cell_area < quantile(seurat\)cell_area, 0.99): Compares each cell’s area to the 99th percentile threshold. The result is a logical vector where each element is TRUE if the corresponding cell’s area is less than the 99th percentile and FALSE otherwise.

seurat[[“SIZE_FILTER_LARGE”]]: Creates a new metadata field named SIZE_FILTER_LARGE in the Seurat object, storing the logical vector.

seurat_CRC1[["SIZE_FILTER_LARGE"]] <- seurat_CRC1$cell_area < quantile(seurat_CRC1$cell_area, .99)

Now we can use ImageDimPlot to visualise the cells which have been flagged for removal.

We can see that these are mostly in the submucosa region.

How do different thresholds behave? Is there a more appropriate one to use? Is any necessary at all?

ImageDimPlot(seurat_CRC1, group.by="SIZE_FILTER_LARGE")

We can use the same approach to create a filter for segmented cells which are very small and likely segmentation arfetacts.

quantile(seurat$cell_area, 0.01): Calculates the 1st percentile of the cell_area values in the Seurat object. This value serves as a threshold to identify the smallest 1% of cells.

seurat\(cell_area > quantile(seurat\)cell_area, 0.01): Compares each cell’s area to the 1st percentile threshold. The result is a logical vector where each element is TRUE if the corresponding cell’s area is greater than the 1st percentile and FALSE otherwise.

seurat[[“SIZE_FILTER_SMALL”]]: Creates a new metadata field named SIZE_FILTER_SMALL in the Seurat object, storing the logical vector.

seurat_CRC1[["SIZE_FILTER_SMALL"]] <- seurat_CRC1$cell_area > quantile(seurat$cell_area, .01)

Now we can use ImageDimPlot to visualise the cells which have been flagged for removal.

We can see that these are more scattered throughout the tissue - but there may be more in the follicular regions.

How do different thresholds behave? Is there a more appropriate one to use? Is any necessary at all?

ImageDimPlot(seurat_CRC1, group.by="SIZE_FILTER_SMALL")

We can check how these values correlate with gene detection rate.

If we filter out small cells, we will remove cells with low numbers of genes detected.

If we filter out large cells, this is not that biased towards overly large counts, as we saw before.

p1 <- VlnPlot(seurat_CRC1, "nFeature_XENIUM", group.by = "SIZE_FILTER_SMALL", pt.size = .1, alpha = .5) + labs(title="Small Cell Filter")

Warning: Default search for "data" layer in "XENIUM" assay yielded no results; utilizing "counts" layer instead.

p2 <- VlnPlot(seurat_CRC1, "nFeature_XENIUM", group.by = "SIZE_FILTER_LARGE", pt.size = .1, alpha = .5)+ labs(title="Large Cell Filter")

Warning: Default search for "data" layer in "XENIUM" assay yielded no results; utilizing "counts" layer instead.

p1 + p2

Adjusting the threshold for what is considered a “small cell” can have significant implications for your analysis, especially in areas with specific cell types such as T-cells, which are small and densely packed in follicular regions. This example demonstrates how changing the threshold to the 10th percentile affects the filtering. In this case, we would probably filter out a lot of good cells that we don’t want to lose! So, be careful when looking at these types of QC metrics!

seurat_CRC1[["SIZE_FILTER_SMALL"]] <- seurat_CRC1$cell_area > quantile(seurat_CRC1$cell_area, .1)

ImageDimPlot(seurat_CRC1, group.by="SIZE_FILTER_SMALL")

Lets set this back to the original 1% threshold.

seurat_CRC1[["SIZE_FILTER_SMALL"]] <- seurat_CRC1$cell_area > quantile(seurat_CRC1$cell_area, .01)

The most important filter is the overall transcript detection. Empty cells or cells with very low transcript count cannot be taken forward for clustering analysis and it is extremely difficult to identify what they may be. Here, we set a threshold of minimum 15 transcripts. This seems quite low - for data from in situ platforms with low noise (Xenium, Merfish, Merscope), this is generally enough to cluster and identify cell types. If your data has more noise (e.g. CosMx), a higher threshold is more appropriate.

seurat\(nCount_XENIUM >= 15*: Compares each cell's transcript count to the threshold of 15. The result is a logical vector where each element is TRUE if the corresponding cell has at least 15 transcripts and FALSE otherwise. *seurat\)TRANSCRIPT_FILTER: Creates a new metadata field named TRANSCRIPT_FILTER in the Seurat object, storing the logical vector.

seurat_CRC1$TRANSCRIPT_FILTER <- seurat_CRC1$nCount_XENIUM >= 15

And we can visualise the cells that we would lose.

We see that we disproportionately would filter out more cells from some regions than others. As pointed out previously, this is likely due to a combination of gene panel coverage in some regions and very small cells in densely packed regions like follicles.

ImageDimPlot(seurat_CRC1, group.by="TRANSCRIPT_FILTER")

Finally, visualizing the counts of negative control codewords, negative control probes, and unassigned codewords helps identify and understand technical artifacts and background noise in your spatial transcriptomics data.

Here, we can see that all control probes and codewords produce yield very little signal, suggesting our data is good quality!

In some cases, high amount of autoflourescence is the cells/tissue can sometimes generate false positive signal and this should be filtered out.

ImageFeaturePlot(seurat_CRC1, "nCount_Negative.Control.Codeword") + scale_fill_viridis_c()

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

ImageFeaturePlot(seurat_CRC1, "nCount_Negative.Control.Probe") + scale_fill_viridis_c()

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

ImageFeaturePlot(seurat_CRC1, "nCount_Unassigned.Codeword") + scale_fill_viridis_c()

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

Although the negative control signal is low, we can nonetheless create a filter to remove cells which have any, although in this case it is probably unnecessary.

seurat_CRC1$PROBE_FILTER <- seurat_CRC1$nCount_Unassigned.Codeword == 0 &
                       seurat_CRC1$nCount_Negative.Control.Codeword == 0 &
                       seurat_CRC1$nCount_Negative.Control.Probe == 0

ImageDimPlot(seurat_CRC1, group.by="PROBE_FILTER")

Finally, we can subset the seurat object based on any/all of the filters we have created earlier.

By combining probe, size, and transcript filters, you can retain only the cells that meet all quality criteria, reducing the impact of technical artifacts and noise on your analysis.

seurat_CRC1 <- subset(seurat_CRC1, PROBE_FILTER & SIZE_FILTER_LARGE & SIZE_FILTER_SMALL & TRANSCRIPT_FILTER)

Warning: Not validating Centroids objectsWarning: Not validating Centroids objectsWarning: Not validating FOV objectsWarning: Not validating Centroids objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating Seurat objectsWarning: Not validating Centroids objectsWarning: Not validating Centroids objectsWarning: Not validating FOV objectsWarning: Not validating Centroids objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating FOV objectsWarning: Not validating Seurat objects

Lets examine the cleaned up object - we have lost a few thousand cells from the analysis.

seurat_CRC1

An object of class Seurat 
541 features across 76132 samples within 4 assays 
Active assay: XENIUM (325 features, 0 variable features)
 1 layer present: counts
 3 other assays present: Negative.Control.Codeword, Negative.Control.Probe, Unassigned.Codeword
 2 spatial fields of view present: COLO_CAN CRC1

Data Normalisation

The SCTransform function in Seurat is used for normalizing single-cell RNA-seq and spatial transcriptomics data. This method models the gene expression counts using a regularized negative binomial regression and removes technical noise while preserving biological variability. The clip.range parameter is used to limit the range of the transformed values, which can help stabilize downstream analyses by limiting the influence of extreme values.

seurat_CRC1 <- SCTransform(seurat_CRC1, assay = "XENIUM", clip.range = c(-10, 10))

Running SCTransform on assay: XENIUM
Running SCTransform on layer: counts
vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
Variance stabilizing transformation of count matrix of size 325 by 76132
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 321 genes, 5000 cells
Found 60 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 325 genes
Computing corrected count matrix for 325 genes
Calculating gene attributes
Wall clock passed: Time difference of 19.44953 secs
Determine variable features
Centering data matrix

  |                                                                                                
  |                                                                                          |   0%
  |                                                                                                
  |==========================================================================================| 100%
Getting residuals for block 1(of 16) for counts dataset
Getting residuals for block 2(of 16) for counts dataset
Getting residuals for block 3(of 16) for counts dataset
Getting residuals for block 4(of 16) for counts dataset
Getting residuals for block 5(of 16) for counts dataset
Getting residuals for block 6(of 16) for counts dataset
Getting residuals for block 7(of 16) for counts dataset
Getting residuals for block 8(of 16) for counts dataset
Getting residuals for block 9(of 16) for counts dataset
Getting residuals for block 10(of 16) for counts dataset
Getting residuals for block 11(of 16) for counts dataset
Getting residuals for block 12(of 16) for counts dataset
Getting residuals for block 13(of 16) for counts dataset
Getting residuals for block 14(of 16) for counts dataset
Getting residuals for block 15(of 16) for counts dataset
Getting residuals for block 16(of 16) for counts dataset
Centering data matrix

  |                                                                                                
  |                                                                                          |   0%
  |                                                                                                
  |==========================================================================================| 100%
Finished calculating residuals for counts
Set default assay to SCT

Principal Component Analysis (PCA) is a dimensionality reduction technique used to identify the primary axes of variation in high-dimensional data. In the context of spatial transcriptomics, PCA helps to reduce the complexity of the data while preserving the most important patterns of variation.

TIP: If your target panel is very small, you can skip this step and carry out clustering analysis directly on gene expression. This can sometimes help with achieving better clustering results.

seurat_CRC1 <- RunPCA(seurat_CRC1)

PC_ 1 
Positive:  CD24, SLC12A2, RRM2, HMGB2, TYMS, PPP1R1B, CA2, CDCA7, SOX9, FERMT1 
       AQP1, STMN1, PCLAF, C1QBP, EPHB3, CMBL, TK1, MKI67, REG4, CEACAM5 
       EGFR, IMPDH2, TUBA1A, SMOC2, S100P, LGR5, CREB3L1, UBE2C, PTTG1, GATA2 
Negative:  IGFBP7, THBS1, TIMP3, DPYSL3, VCAN, IFITM1, CD79A, MAF, CXCR4, CYBB 
       CTSB, ETS1, IL7R, TRAC, TRBC2, CD3E, APOE, ANXA1, CCL5, CD2 
       SERPINA1, SEC11C, RPS4Y1, SOCS3, PLXND1, SPOCK2, MS4A1, RNASE1, CTLA4, CD8A 
PC_ 2 
Positive:  CTSB, APOE, CYBB, RNASE1, SERPINA1, CD14, MS4A7, C1QC, CD163, C1QA 
       IL7R, C1QB, CCL4, FYB1, GPR183, CCL5, CXCR4, CD8A, TRBC2, CD2 
       CD3E, MAF, CD83, TNFSF13B, GZMA, TRAC, CTLA4, PLXND1, CD3D, TIGIT 
Negative:  THBS1, IGFBP7, TIMP3, VCAN, DPYSL3, CDKN2B, CD79A, SEC11C, FZD7, FKBP11 
       DEPP1, RUNX1T1, AQP1, IFITM1, PRDX4, CD24, CKAP4, FRZB, ALDH1B1, TNFRSF17 
       DERL3, CPE, SELENOM, RRM2, SELENOK, TYMS, SLC12A2, PCLAF, EPHB3, SMOC2 
PC_ 3 
Positive:  TRBC2, TRAC, CD3E, CD2, MS4A1, CXCR4, CD8A, CCL5, SPOCK2, CTLA4 
       GZMK, CD3D, GZMA, CD3G, TIGIT, ETS1, KLRB1, IL7R, CD79A, CD6 
       NKG7, CST7, CD5, TRAT1, BANK1, TRBC1, ITK, LTB, CLU, GIMAP7 
Negative:  CTSB, THBS1, RNASE1, IGFBP7, APOE, CYBB, MS4A7, SERPINA1, CD14, TIMP3 
       C1QC, VCAN, CD163, C1QA, PLXND1, C1QB, FZD7, DPYSL3, ANXA1, CDKN2B 
       SOCS3, RUNX1T1, CPE, TUBA1A, FRZB, AQP1, ALDH1B1, TNFSF13B, DEPP1, IL4I1 
PC_ 4 
Positive:  THBS1, CD8A, CCL5, TIMP3, CD2, GZMA, IGFBP7, TRBC2, CD3E, CTLA4 
       TRAC, VCAN, CD3G, CD3D, CCL4, SPOCK2, TIGIT, NKG7, ETS1, KLRB1 
       IL7R, CD6, GZMK, CST7, MAF, AQP1, FOXP3, GNLY, ICOS, CD24 
Negative:  CD79A, SEC11C, FKBP11, MS4A1, TNFRSF17, DERL3, PRDX4, RGS13, SELENOK, CXCR4 
       LRMP, BANK1, CLU, CD79B, SPIB, FCRL1, C2orf88, CYBB, CKAP4, SMIM14 
       GPR183, IRF8, PAX5, CD14, APOE, SELL, TCL1A, MS4A7, RPS4Y1, CD83 
PC_ 5 
Positive:  MS4A1, CLU, BANK1, CXCR4, IGFBP7, SPIB, AQP1, PLXND1, ETS1, DEPP1 
       SOCS3, SMIM14, FCRL1, IFITM1, CPE, IRF8, SELL, ROBO1, LRMP, ARHGAP24 
       IL7R, RPS4Y1, CD83, PAX5, CCR7, TCL1A, LTB, TNFRSF25, MAF, FCER2 
Negative:  SEC11C, CCL5, CD8A, FKBP11, THBS1, GZMA, CD79A, PRDX4, DERL3, TNFRSF17 
       CCL4, NKG7, SELENOK, CD2, CD3E, GNLY, VCAN, CD3G, CKAP4, CD8B 
       CD3D, GZMK, KLRD1, KLRC2, CST7, TRAC, C2orf88, CTSB, CDKN2B, TRBC2 

As before, we can visualise how much variation is captured by each PC.

The ElbowPlot function helps to determine the number of significant PCs to use for downstream analyses. The plot typically shows the amount of variance explained by each PC, and the “elbow” point indicates a natural cutoff.

ElbowPlot(seurat_CRC1, 50)

Plotting the top genes contributing to a specific principal component helps in understanding the biological factors driving the variation captured by that component. This type of plot highlights the genes with the highest loadings, which are the most influential in the principal component analysis.

PC_Plotting(seurat_CRC1, dim_number = 1)

The FeaturePlot function in Seurat is used to visualize the expression of a specific gene across cells in a given dimensionality reduction space (e.g., PCA). This helps to understand how the expression of a gene varies across the principal components.

FeaturePlot(seurat_CRC1, "CEACAM5", reduction = "pca") + scale_color_viridis_c()

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

We can also examine how various PCs are distributed spatially.

Here, we can see that high PC1 loadings enrich in follicular structures and low PC1 loadings enrich in crypt top cells.

ImageFeaturePlot(seurat_CRC1, "PC_1") + scale_fill_viridis_c()

Warning: No FOV associated with assay 'SCT', using global default FOVScale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

We can plot the expression of high (or low) loading genes to visualise how this correlates with our dimensionality reduction.

ImageFeaturePlot(seurat_CRC1, "SMOC2", size=.5) + scale_fill_viridis_c()

Warning: No FOV associated with assay 'SCT', using global default FOVScale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

Next, we will use the reduced dimensionality data for clustering and cluster visualisation.

RunUMAP: Perform Uniform Manifold Approximation and Projection (UMAP) to reduce the dimensionality of the data for visualization. The UMAP plot reduces the high-dimensional data to two dimensions, preserving the local and global structure of the data for visualization. Cells that are close together in the UMAP plot are similar in their gene expression profiles. seurat: The Seurat object. dims = 1:20: Specifies the principal components to use for UMAP.

FindNeighbors: Finding nearest neighbors helps to identify cells that are similar based on their PCA scores, which is used for clustering. seurat: The Seurat object. reduction = “pca”: Specifies that the PCA space should be used for finding neighbors. dims = 1:20: Specifies the principal components to use for identifying neighbors.

FindClusters: Clustering identifies distinct groups of cells with similar gene expression patterns. The resolution parameter controls the granularity of the clustering. seurat: The Seurat object. resolution = 0.7: Sets the resolution parameter for clustering. Higher values lead to more clusters, while lower values lead to fewer clusters.

seurat_CRC1 <- RunUMAP(seurat_CRC1, dims = 1:20)

10:58:18 UMAP embedding parameters a = 0.9922 b = 1.112
Found more than one class "dist" in cache; using the first, from namespace 'spam'
Also defined by ‘BiocGenerics’
10:58:18 Read 76132 rows and found 20 numeric columns
10:58:18 Using Annoy for neighbor search, n_neighbors = 30
Found more than one class "dist" in cache; using the first, from namespace 'spam'
Also defined by ‘BiocGenerics’
10:58:18 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
10:58:32 Writing NN index file to temp file /project/.tmpRsessions/Rtmpyt3o6r/file1d417c7f871a50
10:58:32 Searching Annoy index using 1 thread, search_k = 3000
10:59:24 Annoy recall = 100%
10:59:34 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
10:59:49 Initializing from normalized Laplacian + noise (using RSpectra)
10:59:52 Commencing optimization for 200 epochs, with 3497926 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
11:00:44 Optimization finished

seurat_CRC1 <- FindNeighbors(seurat_CRC1, reduction = "pca", dims = 1:20)

Computing nearest neighbor graph
Computing SNN

seurat_CRC1 <- FindClusters(seurat_CRC1, resolution = 0.2)

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 76132
Number of edges: 2461396

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.9078
Number of communities: 5
Elapsed time: 67 seconds

Next lets visualise the clusters - firstly, based on transcriptome embedding.

DimPlot: Creates a scatter plot of cells in a reduced-dimensional space, by default now using UMAP dimensionality reduction. seurat: The Seurat object containing the dimensionality reduction results and cluster assignments. label = TRUE: Adds cluster labels to the plot. repel = TRUE: Repels the labels to avoid overlapping, making the plot clearer.

DimPlot(seurat_CRC1, label=T, repel=T)

And now lets plot the clusters in tissue space.

We can see that our clusters have quite nice correspondence to distinct spatial regions.

ImageDimPlot(seurat_CRC1, size=.5, axes=T)

Warning: No FOV associated with assay 'SCT', using global default FOV

As before, now we can use Seurat differential expression functions to identify marker genes for specific cell clusters.

FindMarkers: Identifies genes that are differentially expressed in a specified cluster compared to all other cells. seurat: The Seurat object containing the gene expression data and cluster identities. ident.1 = “0”: Specifies the cluster of interest for which marker genes are to be identified. In this case, cluster “0”. max.cells.per.ident = 500: Limits the number of cells to be used from each cluster for the differential expression analysis to 500. This can help to speed up the computation.

markers <- FindMarkers(seurat_CRC1, ident.1="0", max.cells.per.ident=500)

head(markers)

We can visualise expression of cluster specific markers using feature plots

FeaturePlot(seurat_CRC1, "CD3E", label=T, repel=T)+ scale_color_viridis_c(direction=-1)

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

FeaturePlot(seurat_CRC1, "MS4A1", label=T, repel=T)+  scale_color_viridis_c(direction=-1)

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

FeaturePlot(seurat_CRC1, "CEACAM5", label=T, repel=T)+ scale_color_viridis_c(direction=-1)

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

FeaturePlot(seurat_CRC1, "KIT", label=T, repel=T)+ scale_color_viridis_c(direction=-1)

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

FeaturePlot(seurat_CRC1, "CD24", label=T, repel=T)+ scale_color_viridis_c(direction=-1)

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

Or, as in our sequencing ST tutorial, detect and visualise top markers for every cluster.

markers <- FindAllMarkers(seurat_CRC1, max.cells.per.ident = 500)

Calculating cluster 0
Calculating cluster 1
Calculating cluster 2
Calculating cluster 3
Calculating cluster 4

seurat_CRC1

An object of class Seurat 
877 features across 76132 samples within 6 assays 
Active assay: SCT (325 features, 325 variable features)
 3 layers present: counts, data, scale.data
 5 other assays present: XENIUM, Negative.Control.Codeword, Negative.Control.Probe, Unassigned.Codeword, prediction.score.id
 2 dimensional reductions calculated: pca, umap
 2 spatial fields of view present: COLO_CAN CRC1

scCustomize package provides a convenient helper function, Extract_Top_Markers, to extract the top marker genes for each cluster from the output of FindAllMarkers. This function simplifies the process of identifying and retrieving the most significant marker genes for analysis and visualisation.

In this case, we are extracting the top five markers per cluster.

top <- Extract_Top_Markers(markers, num_genes = 5, named_vector = FALSE, make_unique = TRUE)
top

 [1] "SCNN1A"   "LGR5"     "DMBT1"    "EPHB3"    "IL17RB"   "DERL3"    "THBS1"    "TNFRSF17"
 [9] "CPE"      "VCAN"     "MS4A1"    "CTLA4"    "TRBC2"    "CD3D"     "GZMA"     "APOE"    
[17] "RNASE1"   "C1QC"     "C1QB"     "CD163"    "CPA3"     "SLC18A2"  "IL1RL1"   "MS4A2"   
[25] "KIT"     

Clustered_DotPlot function from the scCustomize package provides a convenient and visually appealing way to display expression patterns of top marker genes across clusters using a dot plot. This function not only plots the expression data but also clusters the genes and groups for enhanced visual interpretation. This is an alternative to Seurat DotPlot function.

k = 18: Determines the number of clusters for the hierarchical clustering of genes to enhance visual separation of expression patterns.

We can see that most clusters have unique markers, which suggests the dataset is not over-clustered.

Clustered_DotPlot(seurat_CRC1, features = top, k=18, group.by = "predicted.id")

[[1]]

[[2]]

Additional Spatial Visualisations

The resolution of in situ datasets is typically very high and so it can be difficult to visualise everything in one plot. Below, we will explore different visualisations that can help unpick and understand the data a bit better.

To better visualise spatial distribution of clusters, sometimes it can be useful to subset only certain groups to reduce crowding. Here, we specifically only visualising two selected clusters.

WhichCells: Identifies cells based on specified criteria. seurat: The Seurat object. expression = seurat_clusters %in% c(0, 5): Logical expression to select cells belonging to clusters 0 and 5.

This works with ImageFeaturePlot too. Try it with some genes!

ImageDimPlot(seurat_CRC1, cells=WhichCells(seurat_CRC1, expression = seurat_clusters %in% c(0, 5)))

Warning: No FOV associated with assay 'SCT', using global default FOV

Sometimes, it can be useful to create additional fields of view of the data - for example, zooms of specific regions. First, let’s look at the coordinate system by plotting the data and turning on the plotting of the axes, which are off by default to create nicer looking plots.

This gives us a rough idea on where in the coordinate system to create any subsets or zooms of the data.

For example, if we want to zoom in on the follicle in the top right corner, we can see that it lies roughly between 4000-5000 and 8000-9000 coordinate regions.

ImageDimPlot(seurat_CRC1, axes = T)

Warning: No FOV associated with assay 'SCT', using global default FOV

Cell Type Identification

You can manually annotate your cell clusters, or you can classify them using a reference single-cell dataset. This process is simpler than for Visium data because our data is at the single-cell level, establishing a one-to-one relationship without the need for spot deconvolution.

However, our transcriptome is more limited here, and some cell types may not be well represented. Additionally, our single-cell reference might be missing some cell types that are not well captured by droplet-based technologies but are present in our tissue data.

In this example, we will use a single-cell reference dataset that we prepared earlier.

We will start by reading in the seurat RDS file.

ref <- readRDS("/project/shared/spatial_data_camp/datasets/SINGLE_CELL_REFERENCES/COLON_HC_5K_CELLS.RDS")

Examine the object:

ref

An object of class Seurat 
33556 features across 5725 samples within 3 assays 
Active assay: RNA (33538 features, 2000 variable features)
 3 layers present: counts, data, scale.data
 2 other assays present: HTO, ADT
 2 dimensional reductions calculated: pca, umap

And plot the pre-computed cell clusters. We can see that here we have quite high level annotation.

DimPlot(ref)

We want to evaluate how much structural information is lost in single-cell data when limiting ourselves to the targeted gene set. Accurate cluster prediction is challenging if the current gene set does not adequately identify them. To do this, we will quickly re-embedd the data using only the genes present in our spatial transcriptomics data and keep the original cluster annotations derived from unbiased data.

In this example, we can observe that the limited gene set does a reasonably good job at distinguishing major cell populations. However, it struggles to differentiate between similar cell types, such as myofibroblasts and fibroblasts, as effectively as before.

ref <- SCTransform(ref, residual.features =rownames(seurat))

Running SCTransform on assay: RNA
Computing residuals for the 325 specified features
vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 19008 by 5725
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 5000 cells
Found 80 outliers - those will be ignored in fitting/regularization step

Skip calculation of full residual matrix
Calculating gene attributes
Wall clock passed: Time difference of 28.10518 secs
Determine variable features
Setting min_variance based on median UMI:  0.16
Calculating residuals of type pearson for 319 genes

  |                                                                          
  |                                                                    |   0%
  |                                                                          
  |==================================                                  |  50%
  |                                                                          
  |====================================================================| 100%

Centering data matrix

  |                                                                          
  |                                                                    |   0%
  |                                                                          
  |====================================================================| 100%
Place corrected count matrix in counts slot
Set default assay to SCT

ref <- RunPCA(ref)

PC_ 1 
Positive:  IGFBP7, SOCS3, APOE, IFITM1, TUBA1A, TIMP3, ANXA1, SELENOM, FRZB, VCAN 
       DEPP1, TNFAIP3, THBS1, CXCR4, CCL5, IRF8, RARRES2, CLU, GIMAP7, CD83 
       CD3E, TUBB, GPR183, CCL4, IL7R, PTGER4, ETS1, CD7, RORA, DPYSL3 
Negative:  GUCA2A, SLC26A3, CEACAM7, CA1, CA2, SLC26A2, CEACAM1, AQP8, CEACAM5, CDHR5 
       MS4A12, CA4, SDCBP2, GUCA2B, SELENBP1, MYH14, CD24, MUC12, CD177, CLCA4 
       TMIGD1, HRCT1, CEACAM6, AREG, CES2, C2orf88, GNA11, IL32, SLC6A8, CDKN2B 
PC_ 2 
Positive:  IGFBP7, SOCS3, APOE, GUCA2A, IFITM1, TIMP3, CEACAM7, SLC26A3, TUBA1A, CA4 
       AQP8, CEACAM1, SELENOM, MS4A12, GUCA2B, VCAN, THBS1, DEPP1, CEACAM5, SLC26A2 
       IL32, FRZB, CDHR5, TMIGD1, CLCA4, CD177, SDCBP2, HRCT1, CA1, ANXA1 
Negative:  CXCL3, ADH1C, WFDC2, PPP1R1B, LEFTY1, CCL5, DERL3, CD24, STMN1, UBE2C 
       CD79A, C1QBP, PTTG1, KRTCAP3, CCL20, AREG, TYMS, SLPI, SEC11C, OLFM4 
       CD3E, LGALS2, HMGB2, CDCA7, RHOV, CXCR4, SOX9, GPR183, TKT, LCN2 
PC_ 3 
Positive:  CXCL3, ADH1C, APOE, WFDC2, CD24, SOCS3, IGFBP7, RARRES2, PPP1R1B, LEFTY1 
       SLPI, SELENBP1, KRTCAP3, UBE2C, C1QBP, STMN1, LGALS2, UGT2B17, CA2, LCN2 
       AREG, THBS1, OLFM4, TIMP3, TYMS, NXPE4, CES2, SOX9, CMBL, RHOV 
Negative:  CCL5, CD3E, ANXA1, CXCR4, IL7R, IL32, TNFAIP3, CD7, GPR183, KLRB1 
       CST7, LTB, CD3D, FYB1, CD2, NKG7, LEPROTL1, CCL4, SPOCK2, TRBC2 
       TRBC1, CD8A, CD83, PTGER4, XCL2, CCR7, RORA, CEACAM7, SELENOK, AQP8 
PC_ 4 
Positive:  DERL3, CD79A, TNFRSF17, FKBP11, CCL4, SEC11C, CD83, CA4, AQP8, MS4A1 
       VPREB3, CD79B, GUCA2B, CEACAM7, PRDX4, GUCA2A, CEACAM1, C1QA, C1QC, C1QB 
       DNASE1L3, FCRLA, SELENOM, BANK1, MS4A7, CLCA4, FCER2, LRMP, IL1B, TCL1A 
Negative:  CXCL3, ANXA1, CD3E, ADH1C, CD24, IL7R, WFDC2, CA1, IL32, CA2 
       APOE, CD7, SELENBP1, KLRB1, AREG, PPP1R1B, CCL20, TNFAIP3, RARRES2, SLPI 
       SOCS3, CCL5, LEPROTL1, THBS1, LEFTY1, IFITM1, LGALS2, CD3D, UBE2C, STMN1 
PC_ 5 
Positive:  APOE, CA1, CA2, SLC26A2, SELENBP1, THBS1, SLC26A3, CCL4, VCAN, CD24 
       DERL3, TNFAIP3, IRF8, MUC12, SOCS3, GPR183, CCL5, CD79A, CDHR5, CES2 
       LGALS2, C1QA, TNFRSF17, C1QC, CD14, C1QB, AREG, UGT2B17, ADH1C, FKBP11 
Negative:  CA4, IGFBP7, CLU, RNASE1, GIMAP7, TUBA1A, CA7, BEST4, STMN1, SPIB 
       GUCA2B, WFDC2, OTOP2, HEPACAM2, GUCA2A, AQP1, AQP8, ANXA1, TIMP3, IL32 
       CEACAM6, UBE2C, IFITM1, CEACAM1, KLK1, HES6, RETNLB, TUBB, CLCA4, CTSE 

ref <- RunUMAP(ref, dims=1:20)

14:35:49 UMAP embedding parameters a = 0.9922 b = 1.112
Found more than one class "dist" in cache; using the first, from namespace 'spam'
Also defined by ‘BiocGenerics’
14:35:49 Read 5725 rows and found 20 numeric columns
14:35:49 Using Annoy for neighbor search, n_neighbors = 30
Found more than one class "dist" in cache; using the first, from namespace 'spam'
Also defined by ‘BiocGenerics’
14:35:49 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
14:35:50 Writing NN index file to temp file /project/.tmpRsessions/RtmpONHL0d/file197e2b5e6157d0
14:35:50 Searching Annoy index using 1 thread, search_k = 3000
14:35:52 Annoy recall = 100%
14:36:01 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
14:36:18 Initializing from normalized Laplacian + noise (using RSpectra)
14:36:18 Commencing optimization for 500 epochs, with 231718 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
14:36:35 Optimization finished

DimPlot(ref, label=T, repel=T)

If we visualise the specificity of the gene panel across our single cell reference clusters, we can see that the panel coverage is mainly concentrated across epithelial cells and T-Cells and other immune cells, with few specific markers expressed by stromal cells.

ps <- AggregateExpression(ref, features = rownames(seurat), normalization.method = "LogNormalize", assays="RNA", return.seurat = T)

  |                                                                          
  |                                                                    |   0%
  |                                                                          
  |======                                                              |   9%
  |                                                                          
  |============                                                        |  18%
  |                                                                          
  |===================                                                 |  27%
  |                                                                          
  |=========================                                           |  36%
  |                                                                          
  |===============================                                     |  45%
  |                                                                          
  |=====================================                               |  55%
  |                                                                          
  |===========================================                         |  64%
  |                                                                          
  |=================================================                   |  73%
  |                                                                          
  |========================================================            |  82%
  |                                                                          
  |==============================================================      |  91%
  |                                                                          
  |====================================================================| 100%

ps <- ScaleData(ps, features=rownames(ps))

Centering and scaling data matrix

  |                                                                          
  |                                                                    |   0%
  |                                                                          
  |====================================================================| 100%

pheatmap(LayerData(ps, layer="scale.data"), show_rownames = F)

Next, we can use the standard Seurat integration and cross-classification workflow to transfer single-cell derived labels to our spatial object.

Briefly, the first function identifies anchors between the reference single-cell dataset (ref) and the query spatial dataset (seurat). Anchors are pairs of cells that are considered similar between the datasets. The normalization.method = “SCT” specifies that SCTransform normalization should be used.

The second step transfers the cell type labels from the reference dataset to the query dataset. The anchorset argument specifies the anchors found in the previous step. The refdata = ref$CellType argument specifies the cell type labels from the reference dataset to be transferred. The prediction.assay = TRUE argument indicates that the transferred labels should be stored in a new assay in the query dataset. The weight.reduction = seurat[[“pca”]] argument specifies the dimensionality reduction to be used for weighting the transfer, and dims = 1:30 specifies the number of dimensions to use.

anchors <- FindTransferAnchors(reference = ref, 
                               query = seurat_CRC1, 
                               normalization.method = "SCT")

Normalizing query using reference SCT model
Warning: No layers found matching search pattern providedPerforming PCA on the provided reference using 319 features as input.
Projecting cell embeddings
Finding neighborhoods
Finding anchors
    Found 4981 anchors

seurat_CRC1 <- TransferData(anchorset = anchors, 
                       refdata = ref$CellType, 
                       prediction.assay = TRUE,
                       weight.reduction = seurat_CRC1[["pca"]], 
                       query = seurat_CRC1, 
                       dims=1:30)

Finding integration vectors
Finding integration vector weights
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Predicting cell labels
Warning: Layer counts isn't present in the assay object; returning NULL

Unfortunately, the predicted labels and spatial clusters do not correspond clearly in all cases. This discrepancy is particularly evident in the middle regions of the UMAP, where many cells are predicted as epithelial cells - probably incorrectly!

How to improve this?

Ensure Good Representation of Cell Type Markers in in situ Target Panel Most critically, before undertaking any experiments you want to ensure that there is good representation of all cell types in your target panel - in this case, there is not much to be done as the data has already been generated.

Review and Refine Reference Data: Ensure that the reference single-cell dataset is comprehensive and accurately annotated. If certain cell types are not well represented or annotated in the reference dataset, it can lead to misclassification.

Increase the Number of Dimensions: Increasing the number of dimensions used in the UMAP and PCA steps might capture more variance in the data, leading to better label transfer.

Filter and Preprocess Data: Filtering out low-quality cells or genes and performing additional preprocessing steps can enhance the accuracy of the transfer anchors and, consequently, the label predictions.

Manually Annotate or Correct Predictions: In cases where automatic label transfer is insufficient, consider manually annotating or correcting the predictions for critical regions to ensure accuracy.

As before, we can also visualise the predicted cell labels in tissue space.

ImageDimPlot(seurat_CRC1, cols = cell_colours, group.by = "predicted.id")

Warning: No FOV associated with assay 'SCT', using global default FOV

In line with non-specific predictions, we can also see that the prediction score across these areas is lower.

Outside of stromal cells, we can also see that prediction probability can be low in cells that embedd “between” clusters, for example between core T-Cells and B-Cells, two populations that should be distinct.

This is often the case where cell segmentation is imperfect and partitions transcripts in such a way that it generates “artificial” doublets by pulling in transcripts from an adjacent cell.

FeaturePlot(seurat_CRC1, "predicted.id.score")

For example, if we visualise the lineage markers for T-Cells and B-Cells, we can see that they are often “co-expressed” in the same cells when biologically, they should not be.

The FeatureScatter function in Seurat is used to create a scatter plot showing the relationship between the expression levels of two genes across all cells. This visualization helps to identify potential correlations or patterns between the two genes.

#for the sake of the tutorial we skipped segmentaiton, moved streight to nearest neighbour.

To improve these artefacts, we can try alternative cell segmentation algorithms. What works best is very tissue dependant and there’s no easy one stop solution to this. Cell segmentation algorithms can be divided into a few groups.

Nuclei-based Segmentation algorithms primarily focus on identifying cell nuclei, which are usually more distinct and easier to detect than the cell boundaries. Once the nuclei are identified, the cell boundaries are inferred by expanding around the nuclei. This approach works well in tissues where the nuclei are clearly visible and distinct and in early versions of many in situ platforms, were the only available methods due to only using DAPI stain.

Cell Boundary-Based Segmentation algorithms (e.g. Cellpose) directly segments cells by identifying their boundaries. It is particularly effective for images with complex cell shapes and varying sizes, but this required good cell boundary staining - this is not available for our test dataset. Often cell boundary staining can be non-uniform across different tissues, adding further difficulties. Cellpose version 3 incorporates user-guided model training, which can be very useful for difficult to segment cell types - but this requires time investment to annotate training examples.

Transcript-Density Based Segmentation algorithms, like Baysor segments cells based on the spatial distribution of transcripts. It uses Bayesian inference to assign transcripts to cells, considering both the density and distribution of RNA molecules. This can be very useful for improving cell segmentation where cell boundary stain is not available or not working well.

In this case, we will try re-segmenting our data with Baysor. Here’s the run we prepared earlier - see supplementary material on how to process the data yourself.

baysor <- "/project/shared/spatial_data_camp/datasets/PRECOMPUTED/baysor"

The key output of baysor is the file with transcripts, which have been re-assigned to a new cell identifier.

seg <- read_csv(file.path(baysor, "segmentation.csv"))

head(seg)

#

There will be some transcripts that cannot be assigned to a cell - about 10% in this case. This information is stored under "is_noise" flag. 
This is fairly normal levels of noise.

<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxudGFibGUoc2VnJGlzX25vaXNlKVxuYGBgIn0= -->

```r
table(seg$is_noise)

qplot(seg$assignment_confidence)
table(seg$assignment_confidence > .9)

And transcript confidence - the confidence that the molecule itself is real and not noise.

qplot(seg$confidence)
table(seg$confidence > .9)

We can filter out low confidence and low assignment confidence transcripts here from further analysis. How stringent you want to be depends on whether you want to keep as much data as possible and accept some inaccuracies, or end up with the cleanest possible dataset.

Here, we will filter out transcripts that have not been assigned to cells, and below 0.9 confidence and assignment confidence.

Then, we tabulate a cell by gene matrix from these data.

filtered <- seg[seg$confidence > .9 & seg$assignment_confidence > .9 & !seg$is_noise, ]
mat <- table(filtered$gene, filtered$cell)
mat <- matrix(mat, ncol = ncol(mat), dimnames = dimnames(mat))

Baysor further provides diagnostic info about cells in “segmentation_cell_stats.csv” file, which we will also read in here. The following parameters can be used to filter low-quality cells:

area: area of the convex hull around the cell molecules avg_confidence: average confidence of the cell molecules density: the number of molecules in a cell divided by the cell area elongation: ratio of the two eigenvalues of the cell covariance matrix n_transcripts: number of molecules per cell avg_assignment_confidence: average assignment confidence per cell. Cells with low avg_assignment_confidence have a much higher chance of being an artifact. max_cluster_frac (only if n-clusters > 1): fraction of the molecules coming from the most popular cluster. Cells with low max_cluster_frac are often doublets. lifespan: number of iterations the given component exists. The maximal lifespan is clipped proportionally to the total number of iterations. Components with a short lifespan likely correspond to noise.

stats <- read_csv(file.path(baysor, "segmentation_cell_stats.csv"))
stats <- as.data.frame(stats)
rownames(stats) <- stats$cell
head(stats)

Now, we can assemble a seurat object as before - we first construct a basic object with cell by gene matrix and cell meta data.

seurat_reseg <- CreateSeuratObject(counts = mat, assay = "XENIUM", meta.data = as.data.frame(stats))

Different segmentation algorithms output cell boundaries in various cell formats.

GeoJSON is a format for encoding various geographic data structures using JavaScript Object Notation (JSON). It is widely used for representing spatial features and their attributes and is used by some algorithms to store and output cell segmentation boundaries.

In R, we can read in polygon data from a GeoJSON file using the FROM_GeoJson function.

NOTE: baysor is a 3D cell segmentation algorithm. This means it considered z-stack information. Some algorithms only perform cell segmentation on a representative layer - e.g. Merscope Cellpose takes the middle z-stack. This resegmented data represents 3D segmentations of cells. Since these are 3D segmentations, they may look unusual when visualized as 2D projections.

polygons <- FROM_GeoJson(file.path(baysor, "segmentation_polygons.json"))

In the below code, we extract the polygon coordinates from the data and reformat them into a data frame that Seurat requires to construct a Segmentation object.

polygons <- lapply(1:length(polygons$geometries), FUN=function(x){
  df <- as.data.frame(polygons$geometries[[x]]$coordinates)
  df$cell_id <- paste0("CRef9694c57-", x)
  df
  })

polygons <- do.call(rbind, polygons)
colnames(polygons) <- c("x", "y", "cell_id")
polygons <- polygons[polygons$cell_id %in% Cells(seurat_reseg), ]
polygons <- CreateSegmentation(polygons)

Then, as before, we add both the cell centroid and cell boundaries as segmentations to the seurat object. We skip adding individual molecule coordinates for now.

cents <- CreateCentroids(stats[Cells(seurat_reseg), c("x", "y")])
cents@cells <- Cells(seurat_reseg)
coords <- CreateFOV(coords =list(centroids = cents, segmentation=polygons) ,
                    type = c("centroids", "segmentation"), 
                    molecules = NULL,
                    assay = "XENIUM")

seurat_reseg[["COLON"]] <- coords

From here, we can use the seurat object to visualise various cell meta data - for example, average transcript assignment confidence per cell.

ImageFeaturePlot(seurat_reseg, "avg_assignment_confidence" ) + scale_fill_viridis_c()

Lets filter out low count cells and re-cluster the data as before

seurat_reseg$FILT <- seurat_reseg$nCount_XENIUM >= 15
seurat_reseg <- subset(seurat_reseg, FILT)
seurat_reseg <- SCTransform(seurat_reseg, assay = "XENIUM", clip.range = c(-10, 10))
seurat_reseg <- RunPCA(seurat_reseg)
seurat_reseg <- RunUMAP(seurat_reseg, dims = 1:20)
seurat_reseg <- FindNeighbors(seurat_reseg, reduction = "pca", dims = 1:20)
seurat_reseg <- FindClusters(seurat_reseg, resolution = 0.3)

Visualising clusters, we can see that we already obtain a better separation in the UMAP embedding than before. Though of course, distances in the UMAP space can be very misleading and careful interpretation is required.

DimPlot(seurat_reseg, label=T, repel = T)

Next we visualise the clusters in tissue space.

ImageDimPlot(seurat_reseg)

As before, lets cross-classify our cells using the reference single cell dataset

anchors <- FindTransferAnchors(reference = ref, query = seurat_reseg, normalization.method = "SCT")

seurat_reseg <- TransferData(anchorset = anchors, refdata = ref$CellType, prediction.assay = TRUE,
    weight.reduction = seurat_reseg[["pca"]], query = seurat_reseg, dims=1:30)

Visualising the predictions, we’ve separated T-Cells from B-Cells much better. The stromal clusters still predict poorly, but that is due to poor probe coverage.

DimPlot(seurat_reseg, group.by = "predicted.id")

We can check the distribution in tissue space:

ImageDimPlot(seurat_reseg, group.by = "predicted.id")

FeaturePlot(seurat_reseg,"predicted.id.score")

#skipped segmentaiton, moved streight to nearest neighbour.

Spatial Neighbourhood Analyis

Spatial neighbourhood analysis identifies cells that are spatially close to each other within a tissue section. This technique helps to understand the spatial organization and potential interactions between cells. The same principles used in Visium data can be applied to in situ data.

GetTissueCoordinates: Retrieves the spatial coordinates of the centroids from the Seurat object. which = “centroids”: Specifies that the centroids’ coordinates should be retrieved. rownames(coords) <- coords$cell: Sets the row names of the coords data frame to the cell IDs.

FindNeighbors: Identifies the nearest neighbours for each cell based on their spatial coordinates. as.matrix(coords[, c(“x”, “y”)]): Converts the x and y coordinates to a matrix format. k.param = 20: Specifies the number of nearest neighbours to identify for each cell. return.neighbor = TRUE: Ensures that the function returns the neighbour indices and distances.

TIP: This approach identifies spatial neighbours. If you analysis requires precise identification of directly adjacent or interacting cell neighbours, then a delaunay network based approach would be more appropriate. R package Giotto implements some nice functionalities based on this

coords <- GetTissueCoordinates(seurat_CRC1, which = "centroids")
rownames(coords) <- coords$cell
neighbours <- FindNeighbors(as.matrix(coords[, c("x", "y")]), k.param = 20, return.neighbor=TRUE)

As in our Visium example, we can use it automatically select cells that are adjacent or physically close to some feature of interest. In this case, we want to automatically select any cells that are nearby the group of cells which have been designated as cluster 8, which are mid-crypt epithelial cells.

This can be useful to explore how cells in adjacent cells might influence or interact with each other.

WhichCells identifies the cells that belong to a specific cluster or meet a particular expression criterion

TopNeighbors finds the top n nearest neighbours for a given set of cells based on the k-NN graph. Increasing the n here will return more and more distal spots - tweak to your requirements!

How would you analyse these groups further? What are the adjacent cells? Are they heterogenous? What do they express?

cells <- WhichCells(seurat_CRC1, expression= SCT_snn_res.0.7 == 8)
adjacent <- TopNeighbors(neighbours, cells, n = 10)

Idents(seurat_CRC1) <- "Other Cells"
seurat_CRC1 <- SetIdent(seurat_CRC1, cells = adjacent, "Adjacent Cells")
seurat_CRC1 <- SetIdent(seurat_CRC1, cells = cells, "Cells of Interest")

ImageDimPlot(seurat_CRC1)

seurat_CRC1[["group1"]] <- Idents(seurat_CRC1)

Finding Spatially Correlated Genes

We can use the spatial neighbourhood graph to identify spatially correlated features.

We start by calculating, for each cell, the pseudobulk expression of all cells in it’s local neighbourhood. First of all, lets expand the k.param to include more neighbours - this will control whether we’re including more distal gene expression.

FindNeighbors: Identifies the nearest neighbours for each cell based on their spatial coordinates. as.matrix(coords[, c(“x”, “y”)]): Converts the x and y coordinates to a matrix format. k.param = 50: Specifies the number of nearest neighbours to identify for each cell.

LayerData: Extracts the gene expression data from the specified layer and assay in the Seurat object. layer = “counts”: Specifies that the counts layer should be extracted. assay = “XENIUM”: Specifies the assay from which to extract the data.

as.matrix(neighbours$nn %*% t(mt)): Multiplies the neighbour matrix with the transpose of the gene expression matrix to aggregate the expression levels of neighbouring cells.

neighbours <- FindNeighbors(as.matrix(coords[, c("x", "y")]), k.param = 50)
mt <- LayerData(seurat_CRC1, layer = "counts", assay = "XENIUM")
sum_mtx <- as.matrix(neighbours$nn %*% t(mt))

We can store the neighbourhood-aggregated values in our Seurat object as a separate assay, which we will call “NEIGHBOURHOOD50”. We then normalise the matrix.

seurat_CRC1[["NEIGHBOURHOOD50"]] <- CreateAssayObject(t(sum_mtx))
seurat_CRC1 <- NormalizeData(seurat_CRC1, assay = "NEIGHBOURHOOD50")

We can then apply quick correlation calculations to identify spatially correlated features.

corrgenes <- cor(as.matrix(t(LayerData(seurat_CRC1, assay = "NEIGHBOURHOOD50", layer = "data"))))
diag(corrgenes) <- 0
high_corr_genes <- which(rowMaxs(corrgenes) > .7)
diag(corrgenes) <- 1
heatmap <- pheatmap(corrgenes[high_corr_genes, high_corr_genes], border_color = NA)

Here, we use a very quick cutree function to partition the hierachical clustering into 5 groups of co-localising genes.

TIP- for more advanced module detection, you can use the spatial neighbourhood aggregated matrix with tools such as WGCNA (though it has its own meta-cell functionality).

modules <- cutree(heatmap$tree_row, 5)
modules

Lets visualize some of the detected spatially co-localizing genes. For example, module 1 genes - we can see that CD24 and SOX9 are spatially similar, but not necessarily always expressed by the same cells.

ImageFeaturePlot(seurat_CRC1, "CD24") + scale_fill_viridis_c()
ImageFeaturePlot(seurat_CRC1, "ANXA1") + scale_fill_viridis_c()

#CD24 mosty cancer cells 
# ANXA1 produced by macrophages, neutrophils and epithelial cells 

FeaturePlot(seurat_CRC1, c("ANXA1", "CD24"))

Similarly, module 4 genes which are expressed by different cell types, but show strong co-localisation:

ImageFeaturePlot(seurat_CRC1, "APOE") + scale_fill_viridis_c()
ImageFeaturePlot(seurat_CRC1, "MS4A7") + scale_fill_viridis_c()

We can use Seurat’s module score functionality to add the modules to the meta data as an aggregate score for ease of visualisation.

The AddModuleScore function in Seurat is used to calculate and add scores for predefined gene modules to a Seurat object. These module scores represent the average expression levels of a set of genes (modules) for each cell, which can be used to identify specific biological processes or cell states.

seurat: The Seurat object to which the module scores are added. features = split(names(modules), modules): Specifies the gene modules. The split function is used to create a list of gene sets, where each set corresponds to a specific module. assay = “SCT”: Specifies the assay from which to calculate the module scores. Here, “SCT” refers to the SCTransform-normalized data. nbin = 3: Divides the genes into 3 bins based on their average expression levels to account for differences in gene expression distributions. Need to reduce this from defaults as our number of measured genes is much lower than in single cell analysis. name = “MOD”: Specifies the prefix for the module score names. The resulting scores will be named “MOD1”, “MOD2”, etc.

seurat <- AddModuleScore(seurat, features=split(names(modules), modules), assay = "SCT", nbin=3, name = "MOD" )

Visualising module scores - we can see that we have identified a group of genes co-localising at the base of the epithelial crypts (MOD1) and another module of genes co-localising in lymphoid follicles.

ImageFeaturePlot(seurat, "MOD1") + scale_fill_viridis_c()
ImageFeaturePlot(seurat, "MOD4") + scale_fill_viridis_c()

Detecting Cellular Niches

We can use a similar approach to define cellular niches.

Lets repeat the process of defining spatial neighbourhood expression, with the following modifications:

1. Expand the nearest neighbours even further, to consider 100 nearby cells.

2. Exclude the transcriptomes of the cells themselves and only consider the expression of nearby cells.

This approach will enable us to group cells based not on their identity, but on their microenvironment.

neighbours <- FindNeighbors(as.matrix(coords[, c("x", "y")]), k.param = 100)
diag(neighbours$nn) <- 0 # dont count transcriptome of the cell itself, just neighbours
mt <- LayerData(seurat, layer = "counts", assay = "XENIUM")
sum_mtx <- as.matrix(neighbours$nn %*% t(mt))

How is this useful? Well, now you can cluster cells not on their gene expression values, but gene expression values of surrounding cells. This effectively partitions cells not based on their identity, but on their micro-environment! Using this approach, you can identify tissue niches

Alternative approaches - you could count cell types rather than gene expression values, but that requires you to have finalised cell annotation for your dataset, which is not ideal. So, we do unbiased transcriptomics approach.

How would you run this with cell types?

seurat[["NEIGHBOURHOOD100"]] <- CreateAssayObject(t(sum_mtx))
DefaultAssay(seurat) <- "NEIGHBOURHOOD100"
seurat <- NormalizeData(seurat)
seurat <- ScaleData(seurat, features = rownames(seurat))
seurat <- RunPCA(seurat, features = rownames(seurat))
seurat <- FindNeighbors(seurat, reduction = "pca", dims = 1:10)
seurat <- FindClusters(seurat, resolution = 0.1, cluster.name = "Niches")

Lets visualise the detected “niches”. We can see that we have achieved a coarse partioning of the cells into crypt top, mid-crypt and crypt-base regions, as well as segmenting out follicles and sub-mucosal stroma.

How would you tweak the above approach to generate more or less granular niches?

ImageDimPlot(seurat, group.by = "Niches")

We can tabulate our detected niches with predicted cell type labels (or clusters) to visualise enrichment of different cell types across spatial niches.

For example, as could be expected, T-Cells and B-Cells enrich in Niche 2 (follicular).

comp <- table(seurat$Niches, seurat$predicted.id)
pheatmap(comp, scale="row")

comp <- table(seurat$Niches, seurat$SCT_snn_res.0.7)
pheatmap(comp, scale="row")

Detecting Adjacency Dependent Gene Expression Differences

What if we want to know how expression profiles of specific cell types change based on what they’re adjacent to? We can use the exactly same approach, except for each cell neighbourhood, we only aggregate the expression of cell types of interest.

So now for each cell, our matrix is the sum of expression values of nearby specific cell types.

As an example, here we will ask how expression of epithelial cells changes depending on what they are adjacent to. For simplicity, here we will use predicted cell type IDs from single cell reference data.

How would you select a different group of cells?

neighbors <- FindNeighbors(as.matrix(coords[, c("x", "y")]), k.param = 10)
neighbors$nn <- neighbors$nn[coords$cell, coords$cell]
diag(neighbors$nn) <- 0 # dont count transcriptome of the cell itself, just neighbours?
mt <- LayerData(seurat, layer = "counts", assay = "XENIUM")[, coords$cell]
mt[, seurat$predicted.id != "Epithelium"] <- 0 # Zero counts for cells that are not of interest
sum_mtx <- as.matrix(neighbors$nn %*% t(mt))

As before, we store the aggregated matrix as a new assay, except here we further filter out any cells that do not have any epithelial cells nearby.

seurat[["EPI"]] <- CreateAssayObject(t(sum_mtx[rowSums(sum_mtx) > 0, ]))
DefaultAssay(seurat) <- "EPI"
seurat <- NormalizeData(seurat)

So, how does gene expression of epithelial cells change depending on whether they are in proximity to macrophages (cluster 9), or myofibroblasts (cluster 0).

First, lets visualise the spatial distribution of these cells:

Idents(seurat) <- "SCT_snn_res.0.7"
ImageDimPlot(seurat, cells = WhichCells(seurat,  expression=SCT_snn_res.0.7 %in% c( 9, 0)))

We can use FindMarkers function to identify differences between adjacent epithelial cells for these groups, if we run it on the aggregated EPI matrix:

markers <- FindMarkers(seurat, 0, 9) 
markers

Let’s visualise some of the top hits - we can see they separate nicely in epithelial cells:

VlnPlot(seurat, c("CA7",  "SMOC2"), idents = c(0, 9), assay = "EPI", pt.size = .1, alpha = 0.5)

As a sanity check, if we visualise actual, not epi-adjacent expression in myofibroblasts and macrophages, we can see that there’s hardly any expression at all. So, we are picking up epithelial signal

VlnPlot(seurat, c("CA7",  "SMOC2"), idents = c(0, 9), assay = "SCT", pt.size = .1, alpha = 0.5)

And a visual check

DefaultAssay(seurat) <- "SCT"
ImageFeaturePlot(seurat, "CA7") + scale_fill_viridis_c()
ImageFeaturePlot(seurat, "SMOC2") + scale_fill_viridis_c()

Finally, lets save all of our analysis as an RDS object

saveRDS(seurat, file="colon_in_situ.RDS")

sessionInfo()

###here we plot a percentage proportions of each cell type from integrated datasets 

library(Seurat) 
library(dplyr)
library(ggplot2)
 
meta.data <- CRC_merge[[]]
 
# Calculate percentage
percentages <- meta.data %>%
  group_by(orig.ident) %>%
  mutate(total_count = n()) %>%  # Get total counts per orig.ident
  group_by(predicted.id, orig.ident) %>%
  summarise(count = n(),
            percentage = (count / first(total_count)) * 100)  # Calculate percentage
 
# Plot percentages with labels
ggplot(percentages, aes(x = orig.ident, y = percentage, fill = predicted.id)) +
  geom_bar(stat = 'identity') +
  geom_text(aes(label = round(percentage, 1)),  # Add percentage labels rounded to 1 decimal
            position = position_stack(vjust = 0.5), size = 3) +  # Position labels inside bars
  ylab("Percentage") +
  ggtitle("Percentage of Predicted ID by Orig Ident") +
  theme_minimal()